	/*This program computes logit demand and simulates corresponding price changes assuming Bertrand competition with product differentiation*/

	#delimit;
	set matsize 10000;
	clear;
	clear matrix;
	mata: mata clear;
	local nbrands=7;
	set seed 070910432;
	
	*quietly{;
		set mem 250m;
		set more off;
		do "c:\data\restat_12_9\programs\csolve_do\csolve";
		do "c:\data\restat_12_9\programs\csolve_do\callf";
		do "c:\data\restat_12_9\programs\csolve_do\callf_setup";
		do "c:\data\restat_12_9\programs\csolve_do\antilogit_7_8";
		/*get top level elasticity*/
		mata: mata matuse "c:\data\restat_12_9/analysis_data/oil_top_level";
		mata: st_matrix("ols_top_b",ols_b_top[1,1]);
		mata: st_matrix("iv_top_b",iv_b_top[1,1]);
		mata: st_matrix("ols_top_mean",ols_b_top);
		mata: st_matrix("iv_top_mean",iv_b_top);
		mata: st_matrix("ols_top_vc",ols_vc_top);
		mata: st_matrix("iv_top_vc",iv_vc_top);

		
		
		scalar ols_top_b=ols_top_b[1,1];
		mata: ols_top_b=ols_b_top[1,1];
		
		scalar iv_top_b=iv_top_b[1,1];
		mata: iv_top_b=iv_b_top[1,1];
		
		use "c:\data\restat_12_9\analysis_data\hl_oillogitdata";
		qui tab BRAND2, gen(bra);
		qui tab REGION, gen(reg);
		qui tab date, gen(t);
		qui tab MONTH, gen(m);
		*drop if date>488;
		
		
		local reglist "reg2";
		forval x=3/10{;
			local reglist "`reglist' reg`x'";
		};
		
		su AVOL;
		scalar vt=r(N)*r(mean);
		forval x=1/`nbrands'{;
			su AVOL if bra`x'==1;
			scalar vshare`x'=r(N)*r(mean)/vt;
		};

		forval x=1/`nbrands'{;
			su price if bra`x'==1;
			scalar price`x'=r(mean);
		};
		scalar wavgp=vshare1*price1+vshare2*price2+vshare3*price3+vshare4*price4+vshare5*price5+vshare6*price6+vshare7*price7;
		
		forval x=1/`nbrands'{;
			forval y=2/10{;
				gen b`x'r`y'=reg`y'*bra`x';
			};
		};
		
		forval x=1/`nbrands'{;
			forval y=2/12{;
				gen m`x'r`y'=bra`x'*m`y';
			};
		};
		
		forval x=1/`nbrands'{;
			gen trend`x'=bra`x'*date;
		};
		
		
		forval x=1/4{;
			forval y=2/10{;
				local brandxreg "`brandxreg' b`x'r`y'";
			};
		};
		
		forval x=6/7{;
			forval y=2/10{;
				local brandxreg "`brandxreg' b`x'r`y'";
			};
		};
		
		forval x=1/4{;
			forval y=2/12{;
				local brandxmonth "`brandxmonth' m`x'r`y'";
			};
		};
		
		forval x=6/7{;
			forval y=2/12{;
				local brandxmonth "`brandxmonth' m`x'r`y'";
			};
		};
		
		local trends "trend6 trend7";
		forval x=1/4{;
			local trends "`trends' trend`x'";
		};
		
		local logit "price_plprice bra1 bra2 bra3 bra4 bra6 bra7 `reglist' t2 t3 t4 t5 t6 t7 t8 t9 t10 t11 t12 t13 t14 t15 t16 t17 t18 t19 t20 t21 t22 t23 t24 t25 t26 t27 t28";
		local logit2 "price_plprice bra1 bra2 bra3 bra4 bra6 bra7 `brandxreg' `brandxmonth' `trends'";
		sort BRAND REGION date;
		
		
		sort BRAND REGION date;
		by BRAND REGION: gen T=_n;
		sort T REGION BRAND;
		bysort T: gen id=_n;
		tsset id T;
		
		noisily xi: newey2 lshare_plshare `logit'  if BRAND2~="PRIVATE", lag(1) nocons;
		noisily xi: newey2 lshare_plshare `logit2'  if BRAND2~="PRIVATE", lag(1) nocons;
		
		matrix ols_logitB=e(b);
		matrix ols_logitVC=e(V);
		scalar alpha=_b[price_plprice];
		matrix olsb=_b[price_plprice];
		
		scalar alpha=-1;
		
		forval x=1/`nbrands'{;
			forval y=1/`nbrands'{;
				sca olse`x'b`y'=-price`x'*alpha*vshare`x'+ols_top_b*vshare`x'*price`x'/wavgp;
			};
		};
		
		forval y=1/`nbrands'{;
			sca olse`y'b`y'=price`y'*alpha*(1-vshare`y')+ols_top_b*vshare`y'*price`y'/wavgp;
		};
		
		mat olselas=J(7,7,0);
		forval x=1/`nbrands'{;
			forval y=1/`nbrands'{;
				mat olselas[`x',`y']=olse`x'b`y';
			};
		};
		mat li olselas;
		
		
		
		noisily xi: newey2 lshare_plshare bra1 bra2 bra3 bra4 bra6 bra7 `brandxreg' `brandxmonth' `trends'  (price_plprice=iprice) if BRAND2~="PRIVATE", nocons lag(1);
		matrix iv_logitB=e(b);
		matrix iv_logitVC=e(V);
		
		scalar ivalpha=_b[price_plprice];
		matrix ivb=_b[price_plprice];
		
		
		forval x=1/`nbrands'{;
			forval y=1/`nbrands'{;
				sca ive`x'b`y'=-ivalpha*vshare`x'*price`x'+ols_top_b*vshare`x'*price`x'/wavgp;
			};
		};
		
		forval y=1/`nbrands'{;
			sca ive`y'b`y'=ivalpha*price`y'*(1-vshare`y')+ols_top_b*vshare`y'*price`y'/wavgp;
		};
		mat ivelas=J(7,7,0);
		forval x=1/`nbrands'{;
			forval y=1/`nbrands'{;
				mat ivelas[`x',`y']=ive`x'b`y';
			};
		};
		mat li ivelas;
		
 
		#delimit;

		mat vsbar=J(7,10,0);
		mat pbar=J(7,10,0);
		forval r=1/10{;
			sum AVOL if reg`r'==1;
			sca vt`r'=r(N)*r(mean);
		};
	
	
	
		forval r=1/10{;
			forval b=1/7{;
				qui su AVOL if reg`r'==1&bra`b'==1;
				mat vsbar[`b',`r']=r(N)*r(mean)/vt`r';
				qui su price if reg`r'==1&bra`b'==1;
				mat pbar[`b',`r']=r(mean);
			};	
		};		


		
	#delimit cr;
		mata:
			pbar=st_matrix("pbar")
			vsbar=st_matrix("vsbar")
			
			
			pre_m=I(`nbrands')
			
			
			
			post_m=pre_m
			post_m[4,6]=1
			post_m[6,4]=1
			
			olslogit=J(`nbrands',10,0)
			ivlogit=J(`nbrands',10,0)
			ols_solver_flag=J(10,1,0)
			iv_solver_flag=J(10,1,0)
			olsb=st_matrix("olsb")
			ivb=st_matrix("ivb")
			olsb=-olsb
			ivb=-ivb
			ols_e=-ols_top_b
			iv_e=-iv_top_b
		
			
			for(i=1;i<=cols(pbar);i++){
					pinit=pbar[.,i]
					a=csolve(&anti(),pinit,1e-6,400,vsbar[.,i],pbar[.,i],olsb,ols_e,pre_m,post_m)
					csol=a[1..7,.]
					olslogit[.,i]=csol
					ols_solver_flag[i,1]=a[8,1]
			}
		
		olslogit=olslogit'
		ans=(olslogit-pbar'):/pbar'
		/*take only markets for which we could compute a solution*/
		ols_sim=100*10/(colsum(J(10,1,1)-!rowsum(ans)))*mean(ans,1)
		
			for(i=1;i<=cols(pbar);i++){
					pinit=pbar[.,i]
					a=csolve(&anti(),pinit,1e-6,400,vsbar[.,i],pbar[.,i],ivb,ols_e,pre_m,post_m)
					sol=a[1..7,.]
					ivlogit[.,i]=sol
					iv_solver_flag[i,1]=a[8,1]
			}
		ivlogit=ivlogit'
		ans=(ivlogit-pbar'):/pbar'
		iv_sim=100*10/(colsum(J(10,1,1)-!rowsum(ans)))*mean(ans,1)
	end		

	
	

/***Now that we have our point estimates of costs and approximate simulated price changes, we draw from asymptotic approx of distributions
of estimators and recalculate for each draw.  will look at empirical distribution of these calculated costs and changes to do inference.****/	
#delimit;
drop _all;
/*number of bootstrap iterations*/
local its=1000;
drawnorm `logit2', n(`its') mean(ols_logitB) cov(ols_logitVC);
mkmat _all;


drop _all;
drawnorm e_topbs b2-b23, n(`its') mean(ols_top_mean) cov(ols_top_vc);
mkmat e_topbs;
drop _all;




use "c:\data\restat_12_9\analysis_data\hl_oillogitdata";

tab REGION, gen(reg);
tab BRAND2, gen(bra);
forval x=1/10{;
    mat bpbar`x'=J(`nbrands',`its',0);  
    mat bvsbar`x'=J(`nbrands',`its',0);  
};


forval x=1/10{;
    forval w=1/`its'{;
        preserve;
            keep if reg`x'==1;
            bsample;
            forval y=1/`nbrands'{;
                qui sum price if bra`y'==1;
                matrix bpbar`x'[`y',`w']=_result(3);
                qui sum vshare if bra`y'==1;
                matrix bvsbar`x'[`y',`w']=_result(3);
            };
        restore;
    };
};
#delimit cr;
mata:
    bp1=st_matrix("bpbar1")
    bp2=st_matrix("bpbar2")
    bp3=st_matrix("bpbar3")
    bp4=st_matrix("bpbar4")
    bp5=st_matrix("bpbar5")
    bp6=st_matrix("bpbar6")
    bp7=st_matrix("bpbar7")
    bp8=st_matrix("bpbar8")
    bp9=st_matrix("bpbar9")
    bp10=st_matrix("bpbar10")
    

    ppoint=J(1,10,NULL)
    ppoint[1]=&bp1
    ppoint[2]=&bp2
    ppoint[3]=&bp3
    ppoint[4]=&bp4
    ppoint[5]=&bp5
    ppoint[6]=&bp6
    ppoint[7]=&bp7
    ppoint[8]=&bp8
    ppoint[9]=&bp9
    ppoint[10]=&bp10
        

    bvs1=st_matrix("bvsbar1")
    bvs2=st_matrix("bvsbar2")
    bvs3=st_matrix("bvsbar3")
    bvs4=st_matrix("bvsbar4")
    bvs5=st_matrix("bvsbar5")
    bvs6=st_matrix("bvsbar6")
    bvs7=st_matrix("bvsbar7")
    bvs8=st_matrix("bvsbar8")
    bvs9=st_matrix("bvsbar9")
    bvs10=st_matrix("bvsbar10")
    
    vspoint=J(1,10,NULL)
    vspoint[1]=&bvs1
    vspoint[2]=&bvs2
    vspoint[3]=&bvs3
    vspoint[4]=&bvs4
    vspoint[5]=&bvs5
    vspoint[6]=&bvs6
    vspoint[7]=&bvs7
    vspoint[8]=&bvs8
    vspoint[9]=&bvs9
    vspoint[10]=&bvs10
     

    b=st_matrix("price_plprice")
	e=st_matrix("e_topbs")
    bslogit=J(1000,`nbrands',0)
	flag=J(1000,10,0)	
					
	for(z=1;z<=1000;z++){
		bz=-b[z,1]
		ez=-e[z,1]		
        temp=J(10,7,0)

		for(j=1;j<=10;j++){
            	p=(*ppoint[j])[.,z]
	            pinit=(*ppoint[j])[.,z]
	            vsbar=(*vspoint[j])[.,z]
	            a=csolve(&anti(),pinit,1e-6,400,vsbar,p,bz,ez,pre_m,post_m)
      	        flag[z,j]=a[8,1]
				a=a[1..7,.]
            	
				a=(a-p):/p
	            temp[j,.]=a'
	      }
	bslogit[z,.]=100*10/(colsum(J(10,1,1)-!rowsum(temp)))*mean(temp,1)
	}
end
mata: st_matrix("bslogit", bslogit)
drop _all
svmat bslogit
save "c:\data\restat_12_9\bootstrap_results\ols_oillogit", replace




#delimit;
drop _all;
/*number of bootstrap iterations*/
local its=1000;
	local reglist "reg2";
	forval x=3/10{;
		local reglist "`reglist' reg`x'";
	};
local logit "price_plprice bra1 bra2 bra3 bra4 bra6 bra7 `brandxreg' `brandxmonth' `trends'";
drawnorm `logit', n(`its') mean(iv_logitB) cov(iv_logitVC);
mkmat _all;
drop _all;
drawnorm e_topbs b2-b23, n(`its') mean(ols_top_mean) cov(ols_top_vc);
mkmat e_topbs;
drop _all;




#delimit cr;
mata:
    b=st_matrix("price_plprice")
	e=st_matrix("e_topbs")
    bslogit=J(1000,7,0)
	

	for(z=1;z<=1000;z++){
		bz=-b[z,1]
		ez=-e[z,1]		
        temp=J(10,7,0)

		for(j=1;j<=10;j++){
            	p=(*ppoint[j])[.,z]
	            pinit=(*ppoint[j])[.,z]
	            vsbar=(*vspoint[j])[.,z]
	            a=csolve(&anti(),pinit,1e-6,400,vsbar,p,bz,ez,pre_m,post_m)
      	        a=a[1..7,.]
            	a=(a-p):/p
	            temp[j,.]=a'
	      }
	bslogit[z,.]=100*10/(colsum(J(10,1,1)-!rowsum(temp)))*mean(temp,1)
	}
end
mata: st_matrix("bslogit", bslogit)
drop _all
svmat bslogit
save "c:\data\restat_12_9\bootstrap_results\iv_oillogit", replace
